use "gkbk_nejm_replication_data.dta", replace

* SI Table 1: Demographics (unweighted) *
tab AGE4
tab GENDER
tab RACETHNICITY
tab EDUC5
tab INCOME4
tab dem5
tab gop5

* SI Table 2: OLS Regression *
gen masks_agree = masks
recode masks_agree (5=1) (4=2) (2=4) (1=5)
gen wear_masks = masks_agree
recode wear_masks (5=1) (4=1) (3=0) (2=0) (1=0)
svy:reg masks_agree female ages2 ages3 ages4 black latino education income dem5 gop5 vaccinated transmits2 transmits3
outreg2 using sitable2, word dec(2) label  2aster replace

* SI Table 3: Robustness Checks *
svy:reg masks_agree female ages2 ages3 ages4 black latino education income dem5 gop5 vaccinated transmits2 transmits3
outreg2 using sitable3, word dec(2) label  2aster replace
svy:ologit masks_agree female ages2 ages3 ages4 black latino education income dem5 gop5 vaccinated transmits2 transmits3
outreg2 using sitable3, word dec(2) label append 2aster 
svy:logit wear_masks female ages2 ages3 ages4 black latino education income dem5 gop5 vaccinated transmits2 transmits3
outreg2 using sitable3, word dec(2) label append 2aster 

* Demographic differences mentioned in text *
gen blackorlatinx = 0
replace blackorlatinx = 1 if black == 1
replace blackorlatinx = 1 if latino == 1
* 1-2 weeks after 2nd dose *
svy:mean protections4 if white == 1
svy:mean protections4 if blackorlatinx == 1
svy:reg protections4 blackorlatinx if white == 1 |blackorlatinx == 1
* Unsure *
svy:mean protections5 if white == 1
svy:mean protections5 if blackorlatinx == 1
svy:reg protections5 blackorlatinx if white == 1 |blackorlatinx == 1

* SI Figure 1 *
use "C:\Users\dougk\Box\Atkinson Projects\JAMA letters_postvax_life_outreach\nejm revision\replication files\gkbk_nejmfig1.dta", replace

* Data for Figure 1 *
* Top panel *
use "C:\Users\dougk\Box\Atkinson Projects\JAMA letters_postvax_life_outreach\nejm revision\replication files\gkbk_nejm_replication_data.dta", replace
parmby "svy:mean protections1 protections2 protections3 protections4 protections5", norestore

* Bottom panel *
use "C:\Users\dougk\Box\Atkinson Projects\JAMA letters_postvax_life_outreach\nejm revision\replication files\gkbk_nejm_replication_data.dta", replace
parmby "svy:mean COR10_1-COR10_7", norestore

* Create Figure 1 Panel A -- with data above saved in new file *
* Protection Timing and Outreach NO SIDE EFFECTS*
use "gkbk_nejmfig1.dta", clear
eclplot estimate min95 max95 parmseq if parmseq<13, horizontal ti("Panel A") ylabel(0 "{bf:Protection Timing}" 1 "Immediately after 1st dose" 2 "1-2 weeks after 1st dose" 3 "Immediately after 2nd dose" 4 "1-2 weeks after 2nd dose" 5 "Don't know/Unsure" 6 "{bf:Reported Outreach to Vaccinated Adults}" 7 "Second dose required" 8 "Most effective 1-2 weeks after 2nd dose" 9 "Transmission risk to others unknown" 10 "Continue to wear masks" 11 "Continue social distancing" 12 "Continue to avoid crowds", labsize(small) nogrid notick) xlabel(0 "0%" .1 "10%" .2 "20%" .3 "30%" .4 "40%" .5 "50%" .6 "60%" .7 "70%" .8 "80%" .9 "90%" 1 "100%", labsize(vsmall)) ///
xti("") yti("") graphregion(color(white)) eplot(bar) ciopts(blcolor(black)) estopts(barwidth(0.75) bcolor(gray)) 

* Figure 1, Panel B -- average marginal effects *
use "C:\Users\dougk\Box\Atkinson Projects\JAMA letters_postvax_life_outreach\nejm revision\replication files\gkbk_nejm_replication_data.dta", replace
gen masks_agree = masks
recode masks_agree (5=1) (4=2) (2=4) (1=5)
label var transmits2 "Believe cannot transmit after vax"
label var transmits3 "Unsure if can transmit after vax"
label var income "Income"
label var vaccinated "Vaccinated"
gen wear_masks = masks_agree
recode wear_masks (5=1) (4=1) (3=0) (2=0) (1=0)

gen prediction = .
gen se = .
gen counter = _n
label define vars 1 "Female" 2 "Ages: 30-44" 3 "Age: 45-59" 4 "Age: 60+" 5 "Black" 6 "Latinx" 7 "Education" 8 "Income" 9 "Democrat" 10 "Republican" 11 "Vaccinated" 12 "Believe cannot transmit after vax" 13 "Unsure if can transmit after vax"
label values counter vars
replace counter = . if counter > 13
svy:reg masks_agree female ages2 ages3 ages4  black latino education income dem5 gop5 vaccinated transmits2 transmits3 

margins, dydx(female)
matrix B = r(b)
matrix list B
matrix V = r(V)
matrix list V
replace prediction = B[1,1] if counter == 1
replace se = sqrt(V[1,1]) if counter == 1

margins, dydx(ages2)
matrix B = r(b)
matrix list B
matrix V = r(V)
matrix list V
replace prediction = B[1,1] if counter == 2
replace se = sqrt(V[1,1]) if counter == 2

margins, dydx(ages3)
matrix B = r(b)
matrix list B
matrix V = r(V)
matrix list V
replace prediction = B[1,1] if counter == 3
replace se = sqrt(V[1,1]) if counter == 3

margins, dydx(ages4)
matrix B = r(b)
matrix list B
matrix V = r(V)
matrix list V
replace prediction = B[1,1] if counter == 4
replace se = sqrt(V[1,1]) if counter == 4

margins, dydx(black)
matrix B = r(b)
matrix list B
matrix V = r(V)
matrix list V
replace prediction = B[1,1] if counter == 5
replace se = sqrt(V[1,1]) if counter == 5

margins, dydx(latino)
matrix B = r(b)
matrix list B
matrix V = r(V)
matrix list V
replace prediction = B[1,1] if counter == 6
replace se = sqrt(V[1,1]) if counter == 6

margins, dydx(education)
matrix B = r(b)
matrix list B
matrix V = r(V)
matrix list V
replace prediction = B[1,1] if counter == 7
replace se = sqrt(V[1,1]) if counter == 7

margins, dydx(income)
matrix B = r(b)
matrix list B
matrix V = r(V)
matrix list V
replace prediction = B[1,1] if counter == 8
replace se = sqrt(V[1,1]) if counter == 8

margins, dydx(dem5)
matrix B = r(b)
matrix list B
matrix V = r(V)
matrix list V
replace prediction = B[1,1] if counter == 9
replace se = sqrt(V[1,1]) if counter == 9

margins, dydx(gop5)
matrix B = r(b)
matrix list B
matrix V = r(V)
matrix list V
replace prediction = B[1,1] if counter == 10
replace se = sqrt(V[1,1]) if counter == 10

margins, dydx(vaccinated)
matrix B = r(b)
matrix list B
matrix V = r(V)
matrix list V
replace prediction = B[1,1] if counter == 11
replace se = sqrt(V[1,1]) if counter == 11

margins, dydx(transmits2)
matrix B = r(b)
matrix list B
matrix V = r(V)
matrix list V
replace prediction = B[1,1] if counter == 12
replace se = sqrt(V[1,1]) if counter == 12

margins, dydx(transmits3)
matrix B = r(b)
matrix list B
matrix V = r(V)
matrix list V
replace prediction = B[1,1] if counter == 13
replace se = sqrt(V[1,1]) if counter == 13

gen fig1owerdiff = prediction - (1.96*se)
gen fig1upperdiff = prediction + (1.96*se)
eclplot prediction fig1owerdiff fig1upperdiff counter, xlabel(-1 -.8 -.6 -.4 -.2  0 .2 .4) ti("Panel B") xtitle(" " "Change in Support for Mask-Wearing") ytitle("") ylabel(1 2 3 4 5 6 7 8 9 10 11 12 13,angle(horizontal)) xline(0) horizontal xline(0) graphregion(color(white)) ciopts(blcolor(black)) estopts(color(black))

* SI Figure 1: average marginal effects from logit model *
use "C:\Users\dougk\Box\Atkinson Projects\JAMA letters_postvax_life_outreach\nejm revision\replication files\gkbk_nejm_replication_data.dta", replace
gen masks_agree = masks
recode masks_agree (5=1) (4=2) (2=4) (1=5)
label var transmits2 "Believe cannot transmit after vax"
label var transmits3 "Unsure if can transmit after vax"
label var income "Income"
label var vaccinated "Vaccinated"
gen wear_masks = masks_agree
recode wear_masks (5=1) (4=1) (3=0) (2=0) (1=0)

gen prediction = .
gen se = .
gen counter = _n
label define vars 1 "Female" 2 "Ages: 30-44" 3 "Age: 45-59" 4 "Age: 60+" 5 "Black" 6 "Latinx" 7 "Education" 8 "Income" 9 "Democrat" 10 "Republican" 11 "Vaccinated" 12 "Believe cannot transmit after vax" 13 "Unsure if can transmit after vax", replace
label values counter vars
replace counter = . if counter > 13
svy:logit wear_masks female ages2 ages3 ages4  black latino education income dem5 gop5 vaccinated transmits2 transmits3 

margins, dydx(female)
matrix B = r(b)
matrix list B
matrix V = r(V)
matrix list V
replace prediction = B[1,1] if counter == 1
replace se = sqrt(V[1,1]) if counter == 1

margins, dydx(ages2)
matrix B = r(b)
matrix list B
matrix V = r(V)
matrix list V
replace prediction = B[1,1] if counter == 2
replace se = sqrt(V[1,1]) if counter == 2

margins, dydx(ages3)
matrix B = r(b)
matrix list B
matrix V = r(V)
matrix list V
replace prediction = B[1,1] if counter == 3
replace se = sqrt(V[1,1]) if counter == 3

margins, dydx(ages4)
matrix B = r(b)
matrix list B
matrix V = r(V)
matrix list V
replace prediction = B[1,1] if counter == 4
replace se = sqrt(V[1,1]) if counter == 4

margins, dydx(black)
matrix B = r(b)
matrix list B
matrix V = r(V)
matrix list V
replace prediction = B[1,1] if counter == 5
replace se = sqrt(V[1,1]) if counter == 5

margins, dydx(latino)
matrix B = r(b)
matrix list B
matrix V = r(V)
matrix list V
replace prediction = B[1,1] if counter == 6
replace se = sqrt(V[1,1]) if counter == 6

margins, dydx(education)
matrix B = r(b)
matrix list B
matrix V = r(V)
matrix list V
replace prediction = B[1,1] if counter == 7
replace se = sqrt(V[1,1]) if counter == 7

margins, dydx(income)
matrix B = r(b)
matrix list B
matrix V = r(V)
matrix list V
replace prediction = B[1,1] if counter == 8
replace se = sqrt(V[1,1]) if counter == 8

margins, dydx(dem5)
matrix B = r(b)
matrix list B
matrix V = r(V)
matrix list V
replace prediction = B[1,1] if counter == 9
replace se = sqrt(V[1,1]) if counter == 9

margins, dydx(gop5)
matrix B = r(b)
matrix list B
matrix V = r(V)
matrix list V
replace prediction = B[1,1] if counter == 10
replace se = sqrt(V[1,1]) if counter == 10

margins, dydx(vaccinated)
matrix B = r(b)
matrix list B
matrix V = r(V)
matrix list V
replace prediction = B[1,1] if counter == 11
replace se = sqrt(V[1,1]) if counter == 11

margins, dydx(transmits2)
matrix B = r(b)
matrix list B
matrix V = r(V)
matrix list V
replace prediction = B[1,1] if counter == 12
replace se = sqrt(V[1,1]) if counter == 12

margins, dydx(transmits3)
matrix B = r(b)
matrix list B
matrix V = r(V)
matrix list V
replace prediction = B[1,1] if counter == 13
replace se = sqrt(V[1,1]) if counter == 13

gen fig1owerdiff = prediction - (1.96*se)
gen fig1upperdiff = prediction + (1.96*se)
eclplot prediction fig1owerdiff fig1upperdiff counter, xlabel(-.4 -.3 -.2 -.1  0 .1 .2) xtitle(" " "Change in Probability of Support for Mask-Wearing") ytitle("") ylabel(1 2 3 4 5 6 7 8 9 10 11 12 13,angle(horizontal)) xline(0) horizontal xline(0) graphregion(color(white)) ciopts(blcolor(black)) estopts(color(black))

